(** simplest processing of 2d histogram with slices *)

open Histo
open Histo.H1

let get_sqrt_n h =
    sqrt (integral h)
let se = get_sqrt_n
let dummy_error x = 0.1 *. (abs_float x)
let de = dummy_error

let get_mean h = let (_,mean,_) = max_mean_rms h in mean
let get_rms h = let (_,_,rms) = max_mean_rms h in rms 
let get_max h = let (_,max) = Glob.minmax h.arr in max
let get_min h = let (min,_) = Glob.minmax h.arr in min

(** list of (name,(value,error)) *)
let f_list = [
    "integral",(fun h -> let sum = integral h in sum,(sqrt sum));
    "mean",(fun h -> let mean = get_mean h in mean,(get_rms h)/.(se h));
    "rms",(fun h -> let rms = get_rms h in rms,rms/.(se h));
    "max",(fun h -> let max = get_max h in max,de max);
    "min",(fun h -> let min = get_min h in min,de min);
]

let f_h_list = 
    List.map (fun (name,fn) -> 
        let get h = fn (Histo.get_1d h) in
        (name, get)
    ) f_list

open Glob.Bin

let slice_ranges ~proj ~n h = 
    let axis = match proj with X -> h.H2.xa | Y -> h.H2.ya in
    let range = Axis.range axis in
    let min = Axis.min axis in
    List.map (fun i ->
        let i' = float i and width = range /. float n in
        min +. i' *. width, min +. (i'+.1.) *. width
    ) ( 0--(n-1) )
;;

let get_slice_projection ~proj h2 r =
  let p = match proj with X -> Y | Y -> X in
  H2.projection p (
          match proj with 
            X -> H2.sub ~except_last:true ~x:(Some r) ~y:None h2
          | Y -> H2.sub ~except_last:true ~x:None ~y:(Some r) h2
  )
;;

let process_slices proj ~f ~n hist =
    let h = get_2d hist in
    let slices = slice_ranges ~proj ~n h in
    let proj_hists = 
        Array.of_list (List.map (fun (c1,c2) -> 
            0.5*.(c1+.c2), (get_slice_projection ~proj h (c1,c2))) slices)
    in
    let axis = match proj with X -> h.H2.xa | Y -> h.H2.ya in
    let graph = H1.h {n=n;x1=Axis.min axis;x2=Axis.max axis} 
        { name = h.H2.info.name^(letter proj); 
          title = h.H2.info.title^" slices "^(letter proj) }
    in
    let results_errors = Array.map (fun (c,h) ->
        let res,err = f (H1d h) in
        c,res,err
    ) proj_hists in
    Array.iter (fun (x,y,_) -> H1.fill x y graph) results_errors;
    let errors = Array.map (fun (_,_,err) -> err) results_errors in
    H1.set_sq_errors graph (Array.map (fun a -> a**2.) errors);
    H1d graph
;;

let print_efficiency h1 = 
    match h1 with
        H1d h -> (
            let ineff =  h.arr.(0) /. (H1.weights h) in
            Printf.printf "ineff. is b[0] / sum b:\n";
            Printf.printf "ineff=%f then eff=%f\n" ineff (1.-.ineff);
            flush stdout
        )
      | _ -> print_endline "This operation is only supported fo 1d-histograms"
;;
